arXiv:1507.06364vl [physics.optics] 23 Jul 2015 


Accurate near-field calculation in the rigorous 
coupled-wave analysis method 


Martin Weismann^’^, Dominic F G Gallagher^ and 
Nicolae C Panoiu^ 

^ Department of Electronic and Electrical Engineering, University College 
London, Torrington Place, WCIE 7JE, London, UK 
^ Photon Design Ltd, 34 Leopold Street, 0X4 ITW, Oxford, UK 

E-mail: martin. weismann. 12(9ucl .ac.uk 
24 July 2015 


Abstract. The rigorous coupled-wave analysis (RCWA) is one of the most 
successful and widely used methods for modeling periodic optical structures. It 
yields fast convergence of the electromagnetic far-held and has been adapted 
to model various optical devices and wave configurations. In this article, 
we investigate the accuracy with which the electromagnetic near-held can be 
calculated by using RCWA and explain the observed slow convergence and 
numerical artifacts from which it suffers, namely unphysical oscillations at 
material boundaries due to the Gibb’s phenomenon. In order to alleviate these 
shortcomings, we also introduce a mathematical formulation for accurate near- 
held calculation in RCWA, for one- and two-dimensional straight and slanted 
diffraction gratings. This accurate near-held computational approach is tested 
and evaluated for several representative test-structures and conhgurations in order 
to illustrate the advantages provided by the proposed modihed formulation of the 
RCWA. 
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1. Introduction 

The study of the interaction of electromagnetic waves 
with matter has spawned a large variety of methods 
to analytically or numerically solve the Maxwell equa¬ 
tions (MEs). The results obtained from those meth¬ 
ods are invaluable for understanding, validating, pre¬ 
dicting, and guiding experimental efforts and for the 
design process of electromagnetic and optical devices. 
Amongst the numerical methods used in computational 
electromagnetics, one can fundamentally distinguish 
between time-domain methods, which directly incorpo¬ 
rate the transient behavior of electromagnetic waves, 
such as the finite-difference time-domain method [^, 
and frequency-domain methods, like the finite-element 
frequency-domain method [^, which directly deter¬ 
mine time-harmonic solutions to MEs. The former 
methods are general purpose methods and are capable 
of simulating virtually any electromagnetic structure 
comprising metallic or dielectric objects of arbitrary 
size and shape. However, for certain structures and 
applications, other methods can be superior in terms 
of runtime and accuracy. Eor example, the multiple 
scattering method is a series expansion method in 
the frequency domain that is tailored to calculate inter¬ 
action of light and clusters of spherical particles or 
cylindrical rods [^, and as such it is superior to other 
methods when applied to these particular geometries. 

An important field of applications, where several 
specialized numerical algorithms exist , is the study 
of periodic optical structures, including diffraction 
gratings and periodic metamaterials [^. Amongst 
the numerical methods for periodic structures, such 



frequency domain and is based on the decomposition 
of the periodic structure and the pseudo-periodic 
solution of MEs in terms of their Eourier series 
(ES) expansion, hence the periodicity is naturally 
incorporated into the numerical method. RCWA was 
initially developed for modeling one-dimensional (ID) 
diffraction gratings, but with the introduction of fast 
converging Eourier factorization rules for modeling ID 
and two-dimensional (2D) 


17 18 


19 


21 


structures, 


its formulation for isotropic and anisotropic materials 


22 as well as multilayered and oblique structures 
24 , RCWA has evolved to describe arbitrary, 2D- 
periodic structures. The method has been successfully 
applied to model diffraction gratings, diffractive optical 
elements, surface coatings, spectroscopic applications, 
photonic crystals, and periodic metamaterials. It 
should be noted, that in most cases the functionality 
of these applications of periodic structures depends 


on the electromagnetic far-held, i.e. the propagating 
diffraction orders, and it is known that the RCWA 
delivers fast converging and accurate far-held results, 
as reported in many works 14 - 24 


There is, however, a range of novel applications, 
which rely on the optical near-held of a periodic 
structure, especially at their surface, such as surface- 
enhanced Raman spectroscopy 
harmonic generation 


26 28 


25 


and 


, surface second- 
near-held sensing 


29 32 . These applications require a suitably designed 


near-held distribution, usually optimized for maximum 
held enhancement within specihc spatial domains. Al¬ 
though there have been signihcant advances in experi¬ 
mental optical near-held measurement techniques 33 


these techniques are still in their development state 
and not readily available to accurately characterize 
complex photonic nanostructures. These applications 
and experimental shortcomings lead to a critical de¬ 
mand for numerical methods for periodic structures 
that can facilitate an accurate calculation of electro¬ 
magnetic near-helds, and more importantly, the de¬ 
sign of gratings with optimized near-held patterns. 
With very few exceptions 2Qp4 , a thorough investiga¬ 
tion of numerically calculated near-helds in the RCWA 
has been largely neglected during the development of 
the method, as it is often merely considered a post¬ 
processing step. Moreover, additional reasons for the 
scarcity of reports on the convergence and the accuracy 
of the numerically computed electromagnetic near-held 
in RCWA, are the slow near-held convergence and spu¬ 
rious oscillations displayed by these helds 


34 


In this paper we address the issue of inaccurate 
near-held calculations in the RCWA in several ways. 
Eirst, we use some generic cases of dihraction gratings 
to illustrate the slow convergence of RCWA for near- 
held calculations and reveal the reasons for this 
behavior. Based on this analysis and the continuity 
properties of the electromagnetic helds, a new 
formulation of the held evaluation is proposed, which 
yields faster convergence of the electromagnetic near- 
helds and explicitly fulhlls the continuity properties 
of the electromagnetic held at material interfaces. 
This improved formulation of the numerical evaluation 
of the near-helds is then benchmarked against the 
current formulation, with the aim of making RCWA 
an ehective numerical method for modeling modern 
nanophotonic applications that rely on highly accurate 
near-held calculations. 

The remainder of the paper is organized as follows: 
Section will introduce the reader to the problem 
of inaccurate/spurious near-helds in RCWA for ID 
structures and explain the overall strategy for the 
accurate held evaluation. Section [3] will extend the 
ideas gained from the analysis of ID structures to 
arbitrary, straight or obliquely etched 2D gratings and 
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will present the underlying mathematical formalism. 
Then in section computational results for 2D 
gratings will be presented and discussed. Section 
will investigate near- and far-held convergence of the 
modihed method for slanted gratings, before hnal 
conclusions about the capability of the improved 
RCWA for modeling electromagnetic near-helds will be 
drawn in section |6l 



2. Accurate near-field evaluation for 
ID-periodic structures 


Although the last major roadblock that was preclud¬ 
ing RCWA from becoming a highly effective method 
for modeling ID periodic structures has been removed 
when the correct Fourier factorization rules for trans¬ 
verse magnetic (TM) polarization were introduced 


17 18 , a few topics have continued to attract atten¬ 


tion, namely the convergence for slanted ID-periodic 
gratings 24 and the numerical instabilities associated 
to highly-conductive gratings 35 . These improve¬ 


ments and rehnements of the method are concerned 
with the accuracy of the far-held calculation and they 
achieve excellent convergence results for the coefficients 
of and power carried by the dihraction orders in the 
cover and substrate regions of a grating, namely the 
far-held. This picture changes if the near-held is con¬ 
sidered. Thus, in this section we investigate the numer¬ 
ically calculated near-held inside the periodic grating 
region via RCWA, and explain its slow convergence 
and the origin of the observed high spatial frequency 
oscillations. We also propose an improved held evalua¬ 
tion approach, which yields faster converging near-held 
prohles free of numerical artifacts. 

The ideas developed in this paper are applicable to 
any periodic grating structure with sharp boundaries, 
but in order to conduct a comprehensive assessment of 
the accuracy of the near-held calculations achievable 
with the improved RCWA introduced in this work, 
a number of generic test structures (ID- and 2D- 
periodic, straight and slanted) have been chosen, as per 
hgure[2 To widen the spectrum of test conhgurations, 
three diherent materials are considered for the grating: 
silica (Si 02 ) as a dielectric with low index of refraction 
{ug = nsi 02 — 1-45), silicon (Si) as a high refractive 
index dielectric {ug = nsi = 3.4) and a metal, gold 
(An), with index of refraction Ug = uau = 0-97 + 
1.87i [^. In all examples the substrate is Si02 
{us = nsi 02 )' Normal incidence and TM-polarization 
is considered in all of the sections, i.e. the incident 
electric held amplitude is oriented along the x-axis: 

Ei„c(x) = (l,0,0)^£'i„c(x). 


The example considered in this section is depicted 
in hgure [^a). It consists of a binary grating with 
period, Ai = 1pm, height, h = 0.25 pm, and hlling 


Figure 1. Grating structures mounted on a homogeneous 
substrate (refractive index, Us) considered in this study: a) ID 
binary grating, filling factor, p; b) ID binary grating, slanted 
by 6 = 7r/4 w.r.t. the x-axis. c) 2D grating with rectangle- 
semi-circular cross-section, d) 2D cylindrical grating, slanted by 
6 = tt!A w.r.t. the x-axis. The height of all gratings is denoted 
with h, their period is Ai (and A 2 for 2D-periodic gratings), and 
their refractive index is denoted by Ug. The cover medium is 
vacuum with ric = 1. 


factor, p = 0.5. The grating is illuminated by a 
normally incident plane wave with wavelength, A = 
0.51 pm. 

As the hrst step of our investigation, we used 
RCWA to calculate the far-held. In order to quantify 
the convergence of the RCWA method, the relative 
error of the far-held is dehned as: 


eF{N) = 


^\TiN) - T^ef\2 ^ l^(Ar) _ i^re/|2 


( 1 ) 


where denotes the relative transmitted 

(rehected) power corresponding to a discretization with 
2N + 1 complex FS coefficients. The discretization 
parameter, A, represents the number of harmonics 
retained for each dimension. Moreover, is 

a reference value that is considered to be the exact 
solution or a sufficiently good approximation of the 
exact solution. Due to the absence of the exact solution 
of the dihraction grating problem, the reference values 
are chosen to be numerical values obtained by high-A 
simulations; in our case we chose and 

ppref _ /^(905)^ 

The far-held relative error, eF{N)^ for increasing 
number of harmonics, N = 5,..., 640, is depicted in 
hgure As this hgure illustrates, RCWA converges 
quickly for all three materials. Specihcally, in 
order to achieve a self-error of ei?(A) < 1% (as a 
generic criterion adopted here for an accurate far- 
held calculation), for the three materials, silica (red 
crosses), silicon (blue stars), and gold (green triangles), 
a relatively small number of harmonics is necessary, 
namely A > 5, A > 25, and A > 13, respectively. 

As mentioned in the introductory section, the 
far-held of optical gratings and periodic structures 
has been the physical quantity of most interest from 
experimental point of view, hence the characterization 
of a numerical method by means of the far- 
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numerically obtained using N harmonics: 




(3) 


Here, denotes the FS reconstruction given by the 
2N + 1 central Fourier coefficients, Ean^ which are 
calculated by RCWA: 


N 

n=-N 



Figure 2. Far-field relative error vs. the number of harmonics, 
determined for three different materials. 

field convergence has usually been the adopted 
strategy. This approach, however, largely neglects 
the electromagnetic near-held predicted by a specihc 
method. On the other hand, the near-held is 
of fundamental importance for modeling plasmonic 
ehects or optical nonlinear phenomena in devices with 
size comparable to or smaller than the operating 
wavelength, ehects whose description relies on accurate 
calculations of the electromagnetic near-held. 



Figure 3. Self convergence of the tangential and normal electric 
field components Ez and Ex and Ex inside the grating structure 
described by their self-errors AEt{N) (green triangles), AEn{N) 
(red crosses) and AEn{N) (blue circles), respectively for the 
three benchmark structures made of Silica (a), Silicon (b) and 
Gold (c). 

In order to characterize the numerically obtained 
near-helds, we dehne the grating norm, H-H^, of a scalar 
or vectorial function, /, in the grating region as follows: 

fh rAi/2 \ 

/ / \f{x,z)\^ dx dz\ , (2) 

Jo J-A 1 I 2 J 

where the ^-integration extends over the bulk of 
the periodic region. The grating norm is used to 
dehne the near-held error, of the scalar held 

components E^\ a = of a near-held, 


Similarly to the far-held calculations, = 

E(905) obtained by a high-resolution RCWA 
simulation with N = 905. The near-held self¬ 
convergence for the tangential component, Et = 
Ez^ and the normal component, E^ = E^^ of the 
electric held is depicted in hgure For all three 
materials, the self-convergence of Et (blue circles) is 
fast and comparable to the far-held self-convergence 
(see hgure [^. The normal component, E^ (red 
crosses), however, exhibits much slower convergence 
and even at the highest numerical resolution of A" = 
640 a relative error of > 0.9% still remains, 

whereas the error of the tangential held, < 

0.08%, for all materials. 

One intriguing question raised by the data plotted 
in hgures and is why the normal component of the 
near-held converges much more slowly than the far- 
held power and the tangential component of the near- 
held. Or put it the other way around: how can the 
far-held converge quickly when the near-helds have not 
converged yet? There are two factors that explain this 
behavior: i) The far-held consists of a superposition 
of a small number of propagating dihraction orders, 
i.e., plane waves. Hence, the far-held requires a small 
number of FS components to be reconstructed and does 
not suher from Gibb’s phenomenon, ii) RCWA does 
not depend on the representation of the discontinuous 
normal held En. Instead, the method relies on 
the correct Fourier factorization of the continuous 
normal component of the displacement held, which 
can be accurately described by a FS, i.e. without 
spurious oscillations. The second observation is at 
the core of the accurate near-held calculation that 
is introduced in this work. Specihcally, instead of 
reconstructing a discontinuous physical quantity, i.e. 
A^(x), directly from its FS coefficients, it is more 
ehective to reconstruct a continuous quantity, the 
normal component of the displacement held, Dn{'x.)^ 
and then divide it by a discontinuous quantity, the 
electric permittivity, 5o^(x), which is known in the 
space-domain where one seeks to solve the dihraction 
problem. 

For ID-periodic gratings, we dehne the modihed 
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conventional 

N=2l 




-A/2 0 yl/2 


♦ b) 



Figure 4. Near-field distribution of the normal component of 
the electric field, \En\ = \Ex\^ in and around the Au grating. 

a) Calculated with N = 21 harmonics using the conventional 
RCWA method (unphysical field oscillations can be observed). 

b) Calculated using the improved formulation, En = Dnis, 
and A” = 21 harmonics, c) Calculated using the improved 
formulation and N = 640 harmonics. 

normal component of the electric field: 

= £0 ^ A^4x)/e(x), (4) 

where is the vacuum permittivity and ^(x) the 
relative permittivity at position, x. Note that E^\^) 
and En^\'K) represent the same physical quantity, 
namely the normal component of the electric field. 
However, whereas the former is found by using RCWA 
to solve directly for the electric field, the latter one is 
determined by first calculating the displacement field 
and then the electric field via equation @- 

The error, of the modified normal 

component, En^ is shown in figure and encoded in 
green triangle. It is found to converge as fast as the 
fast-convergent tangential component, Et^ and as fast 
as the power in the far-held shown in hgure Even 
at the highest considered resolution, N = 640 , the 
conventional formulation of the RCWA yields a self 
error of > 9 • 10 “^ for all three materials. 

By contrast, the same self error e{En^^} < 9 • 10 “^ of 
the modihed normal held En can be achieved by using 
as few as A" = 70 harmonics. 

The spatial prohle of the electric held in the 
grating region illustrates the full benehts of the 
modihed held calculation. Figure depicts 

the conventional normal held component lEn'^^l = 



Figure 5. Closeup of the normal component of the electric 
field, \En\ and |Tn|, near the metal-air interface at 2 ; = h/2, 
computed using the conventional and modified RCWA methods, 
respectively. 

\Ex \ in the gold grating for a moderately coarse 
discretization of A = 21 harmonics. It can be seen 
that I An \ exhibits unphysical oscillations with a 
spatial frequency equal to the period of the smallest 
spatial frequency component in the FS expansion of the 
solution. This is the well known Gibb’s phenomenon, 
which occurs when describing a discontinuous function 
with a truncated FS. On the other hand, the modified 
normal field. An, does not suffer from such spurious 
oscillations at the interface. In particular, even for a 
small number of harmonics, A = 21, An ^ is smooth, 
as per figure [^b). At very large number of harmonics, 
A = 640, the modified normal field is free of any 
numerical artifacts, as can be seen in figure [^c). 

The improved formulation of RCWA exhibits 
another benefit, namely An is by construction 
discontinuous and exactly fulfills the corresponding 
boundary condition, 

e(in)E(in)(x,) ■ n = (x.) ■ n, (5) 

at surface points of the grating, x^, where and 

^(in) ^(out)^ denote the electric field and 

permittivity inside (outside) the grating, respectively, 
and n is the unit vector normal to the surface. In 
the conventional formulation of the RCWA, the field 
E^^ does not satisfy equation ([^, because E^^ is 
- as a FS containing a finite number of terms - 
inherently continuous. The closeup of the interfacial 
field around x = 0.25A in figure emphasizes 

these ideas. Thus, the normal component of the 
field calculated using the conventional formulation of 
RCWA shows spurious oscillations for both small A = 
21 (dashed green) and high A = 640 (solid red) number 
of harmonics, whereas the modified formulation is free 
of such unphysical oscillations and discontinuous for 
any number of harmonics, A (see the dotted purple 
and dashed-dotted blue lines corresponding to A = 21 
and A = 640, respectively). 
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It should be clear now that the modified normal- 
field evaluation in ID by means of the displacement 
field represents an improvement of the conventional 
evaluation of En in two ways: First, it exhibits optimal 
self-convergence in the sense that it is as accurate as the 
far-field and the continuous tangential field component. 
Second, it explicitly fulfills the boundary condition 
(§ at the material interfaces. Equally important, 
the improved near-held calculation is achieved with 
minimal additional computational cost and does not 
alter the mathematical framework of the core RCWA 
algorithm in ID. 

3. Formulation of accurate near-field 
evaluation for 2D-periodic structures 

The ideas of the previous section are straightforward 
to implement because in ID-periodic structures there 
is a trivial and unambiguous distinction between the 
tangential and normal components of E. However, 
this situation becomes more intricate in the case of 
2D-periodic diffraction gratings. 

As it will be shown now, accurate near-held 
evaluation is strongly related to correct Fourier 
factorization in 2D-periodic structures. Fourier 
factorization means in this context the decomposition 
of a product of periodic functions, / = ^ • h, into its 
periodic factors g and h. Depending on the continuity 
properties of the factors and the product, diherent rules 
must be applied in order to obtain fast convergence 
when increasing the number of FS terms. According 
to these rules, if / and h possess simultaneous 
discontinuities, but g is continuous at those locations, 
the product rule yields fast convergence with respect 
to the number of FS terms, namely one uses f = g • h. 
Moreover, if g and h are simultaneously discontinuous, 
but / is continuous, the inverse rule should be used, 
i.e., / must be factorized as / = (1/^)“^ -h. A rigorous 
explanation of these rules and the solution to the ID 
Fourier factorization problem is given in [^. 

Because of the trivial distinction between the 
continuous (tangential) and discontinuous (normal) 
components of the electric held in ID, Fourier 
factorization is straightforward in that case. But for 
2D-periodicity, three diherent approaches to achieve 
the correct Fourier factorization have been proposed: i) 
Approximate the material boundaries by a coordinates- 
aligned staircase-contour [^. ii) Devise a coordinate 
system in which a given grating geometry is coordinate 
system aligned and use approach i to obtain the 
correct Fourier factorization 37 . Hi) Construct 


a normal vector held (NVF) to decompose E into its 
normal and tangential components and then apply the 


corresponding correct factorization rules to them 21 


the decomposition of the electric held into normal 
and tangential components, it is natural to use the 
factorization approach (iii)^ the NVF approach. It 
is out of the scope of this paper to derive the full 
formulation of 2D-RCWA, so that only the crucial and 
unconventional steps will be given here. Thus, the 
normal and tangential components of D = ^E have 
to be decomposed using the product- and inverse rule, 
respectively. This leads to the following relations for 
the Qf-component of the displacement held 21 24 38 


[Do] = £0 F 

/3=1 


(6) 


where 5o,f3 is the Kronecker delta. Here, [/] denotes 
the vector of FS coefficients of a scalar function 
/, 1^] is the Toeplitz matrix of FS coefficients 
of and the matrix is given by = 

pA[iV„iV^] + [7V„7V^]A) with A = [£] - [1/e]-' 
and is the ^-component of the NVF, N = 
(Vi, V 2 , Vs)^, of the material boundary. The matrix 
~ ^oi(3 implements the three steps of the 
Fourier factorization: the decomposition of E into 
normal and tangential components {Ng in {NcNg}), 
factorization using either the inverse rule (|l/^]~^) 
or the product rule (I^J), and back-projection to 
Cartesian coordinates {N^ in |Vq,V^]). 

Note that one can also use asymptotically equiva¬ 
lent dehnitions of However, our investigations 

have shown that despite the fact that the choices 
Aa /3 = AlA^A/j] and Aap = [AaAT/jJA produce 
similar convergence speed, they do not conserve the 
power for lossless structures, whereas the choice A^g = 
|Vck]A|V/ 3 ] yields power conservation but at the price 
of slower convergence. All three formulations are 
asymptotically equivalent with respect to the number 
of terms in the FS expansion, due to the commutation 


of Toeplitz operators 38 


Normal vector fields can be constructed analyti¬ 
cally for a variety of structures and automated algo¬ 
rithms to obtain a NVF for arbitrary grating geome¬ 
tries have been developed 39 . It should be noted that 


this formulation allows inclined NVFs, i.e. NVFs with 
simultaneously non-vanishing x, and z components. 
Such a NVF is required to accurately model obliquely 
etched structures (see figure Bd)) and hence can be 
viewed as a generalization of the methods presented 
in 1^ {Ny = 0, figure [^b)) and {N^ = 0, fig- 
ure]Tpc)). 

Given the constitutive relation § , one can derive 
the RCWA eigenvalue-problem 


M 


Since the accurate near-held evaluation relies on 


[ 5 .] \ 

[Sy] 

m 

[Uy] j 


[ 5 .] \ 

[Sy] 
[U.] 
[Uy] / 


( 7 ) 
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for the electromagnetic modes described by [S] and 

E and propagation constant, /3. The most general 
mulation of the system-matrix M reads: 



Kx 

Kx S/\zy 

S/S.ZX 

KyBAzy 

^yx 

T ^yz^^zx Kx Ky 

Kx Kx “h AyzBAzy C.y 

\ Lee 

“ AxzBA^x ^y^y 

Kj/Kx Axz^Azy Axy 



KxBK^ l-KxBKx\ 

K^BK^-I -K^Kx 

AyzB\^y AyzB\^x 

AxzBK^ AxzBKx J 


with B = ([e]] - ^ and C„ = |e] - The 

matrices = diag {kan) are diagonal matrices of the 
in-plane propagation constants, k^n for a = x and 
kyn for a = y, of the diffraction orders and I is the 
identity matrix. The size of matrices B, Cq,, Kq,, and I 
is Nq X Nq^ whereas the size of M is 4A^o x 4A^0: where 
No = { 2 N + 1)^ is the total number of FS coefficients 
in the calculation and the same number of FS terms 
for the truncation of the FS in the x- and ^-direction 
is assumed. 

The mode amplitudes [S] and [U] are determined 
by matching the tangential components of the 
electromagnetic field at the top and bottom of 
the grating region, whereas multilayered structures 
are described by the staircase approximation in 
conjunction with the numerically stable 5-matrix 
algorithm. This yields the FS coefficients, [E], of the 
electric field everywhere in the grating region, in the 
cover, and the substrate. 

Let us now denote by 7Z{[f]) the FS reconstruction 
of a Fourier coefficient vector [/], 


where 1 denotes the 3x3 identity matrix and NN^ 
is the 3x3 projection matrix defined by the NVF 
at any point in space, x. Note that by construction, 
[D^] and [E^] are FS coefficient vectors of vector fields 
that are continuous at material interfaces. Hence their 
reconstructions, 7l{[Dn]) and 7^([Et]), do not suffer 
from Gibb’s phenomenon at the interface. With this 
observation in mind, the electric field in the improved 
RCWA method for 2D-periodic structures at a point, 
X, is given by: 

eW(x) =£o'7^([D„])(x)/e(x) +7^([Et])(x) (10) 

and is expected to yield fast near-held convergence and 
non-oscillatory spatial held prohles. To investigate the 
validity of these predictions, the accurate held eval¬ 
uation was implemented in a commercially available 
RCWA computer software, OmniSim/RCWA [io] . 

It should be noted that the other electromagnetic 
helds can be easily calculated with our improved 
method, too. Specihcally, the displacement held, D, 
can be evaluated using the modihed electric held £’, 
namely D^^^(x) = 5o^(x)E^^^(x), and hence will have 
the same convergence properties as E. The magnetic 
induction B and the magnetic held H = /iB do not 
require special attention, because they are continuous 
in non-magnetic materials and hence behave similar 
to the continuous tangential component of the electric 
held. 

4. Quantification the accurate near-held 
evaluation in 2D-periodic structures 


= X] Yi fnmexp(im—x + in—y 
n=—N m= — N ^ ^ ^ 

where fnm — [/](n+Ar)(2Ar+l)+m+Ar+l ^0 FS 

coefficients of /. Within the conventional RCWA 
framework, each component of E is evaluated as: 

EW=7^([E«]). (8) 

This relation does not take into account the continuity 
properties of the diherent components of E and hence 
will lead to spurious oscillations and slow convergence 
of the near-held as it was seen in section [21 

The modihed held evaluation for 2D-periodic 
dihraction gratings requires one to dehne the continu¬ 
ous normal component of D and the tangential com¬ 
ponent of E. Their FS coefficient vectors are given by: 


\ In this section, the improved formulation for accurate 
J ’ near-held calculations in 2D-periodic structures is 
assessed using two test structures under diherent 
conhgurations. To this end, we hrst extend the 
dehnition of the grating norm <§ of a scalar or vector 
function, /, to 2D-periodic structures in the following 
straightforward way: 

/ ph pA2/2 pAi/2 

II/IIg = / / / \f{x,y,z)\‘^ dx dy dz\ , 

\Jo J-A 2 I 2 J-A 1/2 J 

( 11 ) 

where the integral is evaluated over the three- 
dimensional grating region. 

4 . 1 . Analysis of a ID-periodic grating using 
2D-RCWA 


[Dn] = \eo ([NN^lIl/er' + [1/er'lNNq) [E], 

(9a) 

[E*] = [l-NNq[E], (9b) 


The hrst 2D-periodic dihraction grating under consid¬ 
eration is a ID binary grating, as shown in hgurej^a), 
rotated by 7r/4 in the x — ^-plane. It should be obvious 
that it can be modeled as a double-periodic 2D grating 
with periods Ai = A 2 = v^A = ^/2 pm, where A = Ai 
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Figure 6. Computational results for a rotated binary ID 
grating, a) Error of calculated far-field vs. determined 

for three gratings made of different materials. The decrease 
in the error follows that of the far-held (dashed lines) of the 
ID simulations (see also section]^, b), c), d) Near-held error 
corresponding to the silica, silicon, and gold grating, respectively. 


is the period of the grating when it is viewed as a ID- 
periodic structure. For clarity, the primary unit cell 
of the grating is depicted in the inset of figure [^a) . 
With this choice, the reference quantities ^ 

and can be calculated using ID simulations, an 


approach inspired by an example in 21 


The error in the calculation of the far-held, ep 
from equation 0, when 2D simulations are employed 
is depicted in hgure [^a). The convergence of the 
calculations for the silicon and gold gratings closely 
follows the convergence trends observed when ID 
simulations are performed and, as expected, it shows 
somewhat worse, yet still good, agreement for the silica 
structure. This difference is explained by the fact that 
the NVF introduced in the 2D-RCWA formulation is 
discontinuous away from material interfaces and hence 
it can degrade the convergence rate, especially for the 
low-index of refraction contrast case. 

The conclusions of our analysis of the convergence 
of the near-held are summarized in hgures [^b)-[^d) . 
Thus, for both the silica and gold dihraction gratings, 
the normal component of the modihed electric held, 
E, exhibits faster convergence than this same held 
component calculated using the conventional form 
of RCWA. In both cases, e{En^^} is one order of 



Figure 7. a) Normal component of the electric field, \En\: in the 
grating region at 2 : = h/2, determined by using the conventional 
RCWA and N = 27. b) Normal electric field component, \En\, 
determined for the same grating parameters as in a) but using 
the improved algorithm. The blue vertical line was added for 
clarity and merely separates the two plots. 


magnitude smaller than For the silicon 

grating, only marginal diherences between the two 
formulations can be observed. This is in agreement 
with the results obtained in the ID case, as per 
hgure l^b). For small and E*^^^ are 

determined with a comparable degree of accuracy, 
but for the larger number of harmonics considered in 
hgure [^c), i.e. N = 30, a higher accuracy of our 
improved formulation of the RCWA can clearly be 
observed. 

These conclusions are further validated by the 
prohle of the electric held, as presented in hgure 
This hgure shows the spatial distribution of the normal 
component of the electric helds, E^^^^ and E^^^\ 
calculated in the median plane of the grating. A simple 
examination of these held prohles conhrms that the 
spurious oscillations of the held near the surface 

have much smaller amplitude as compared to that of 
the variations of E^^^\ Moreover, a closer inspection 
of the surface-helds shows that the boundary condition 
^ is fulhlled by E^^^^ only. This hrst test-case already 
reveals that the improved near-held evaluation is more 
accurate in the case of 2D-periodic structures, too. 

4 . 2 . Near-field calculations for an intrinsically 
2D-periodic grating 

In order to thoroughly test the near-held evaluation 
for 2D-periodic structures using the improved RCWA 
presented in this article, in what follows we consider 
the challenging test structure depicted in hgure [^c). 
The grating region consists of a coordinate system 
aligned parallelepiped with the length of the sides 
aligned to the x-, ^-, and z-axis being a = 0.5Ai, 2a, 
and h, respectively, placed adjacently to a semicircular 
cylinder with radius a and height h, with A = Ai = 
A 2 = 0.25 pm. The structure is illuminated normally 
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Figure 8. Computational results for the 2D grating shown in 
figure [^c). a) Error of calculated far-field vs. A/”, determined 
for three gratings made of different materials, b), c), d) Near¬ 
field error corresponding to the silica, silicon, and gold grating, 
respectively. 

by a x-polarized plane wave with wavelength A = 
0.5 pm. 

Since it is generally computationally time consum¬ 
ing to obtain high-accuracy solutions in the case of 
2D-periodic structures and due to the fact that the 
higher the ratio A/A and the refractive index |n|, the 
more harmonics are necessary to achieve convergence 
[M] , a relatively small period-to-wavelength ratio of 
A/A = 0.25 pm/0.5 pm = 0.5 is chosen for this exam¬ 
ple. 

As reference values in the definition of the far-held 
error, ep from equation Q, and = 

namely results obtained from simulations with 
N = 31 are chosen. The convergence of the far- 
held, ep^N)^ is shown in hgure [^a), where the far- 
held physical quantity considered are the transmission 
and rehection coefficients, T and i?, respectively. As 
expected, the fastest convergence can be observed for 
the silica grating, because it has a low refractive index. 
The numerical error obtained for the gratings made of 
gold and silicon are one and two orders of magnitude 
larger than in the case of the silica grating, respectively. 
This behavior is similar to that seen in the ID case for 
a small number of harmonics, A" < 30 (cf. hgure [^. 

Figures |^b)-[^d) contain the dependence of the 
near-held self-error determined for the three material 
conhgurations, silica, silicon, and gold, respectively. 


The three physical quantities plotted in each case are 
the ^-component of E and the in-plane components, 
^xy •= -^ 7/5 0)^ and ^xy •= {Ex^Ey^O)'^, 

obtained by using the conventional and improved 
versions of RCWA, respectively. Note that the in¬ 
plane component contains the discontinuous normal 
component, N-E (N-E), and the continuous tangential 
component, (1 — NN^)Ea,^ ((1 — NN^)Ea,^). 

It can be seen that in all cases, the ^-component 
of E, which is continuous at vertical surfaces inside 
the grating region, converges much faster than the 
in-plane component, in both the conventional and 
modihed formulations. In addition, the modihed 
formulation leads to a somewhat smaller error than 
the conventional formulation. Specifically, it was found 
that e{El^^} 0.5e{Ei^^} for the silica grating 

and ^ 0.9e{Ei^^} for the silicon and gold 

gratings. However, the corresponding convergence 
speed, i.e. the slope of and is 

the same. Moreover, in the modified formulation of 
the RCWA the convergence speed of the tangential 
component of the near-held is larger than that of 
the in-plane component. Three factors contribute 
to this behavior: i) The decomposition of the near- 
held in a normal and tangential component in 2D- 
periodic structures relies on the specihc dehnition of 
the normal vector held, N and it is not directly 
performed in Cartesian coordinates as in the ID case. 
Hence, the inexact held decomposition by N introduces 
an additional error, ii) The held of normal vectors 
characterizing the structure is only uniquely dehned 
at the interfaces dehning the grating, except at the 
corners, and, more importantly, away from the grating 
surface. This ambiguity allows and can lead to choices 
of NVFs, which are not optimal for the convergence 
and accurate calculation of the near-held. Hi) The 
normal vector held itself has discontinuities, which can 
cause additional oscillations in the spatial prohle of the 
electromagnetic held. 

It is also worthwhile to investigate the spatial 
prohle of the near-held. The dominant x-component of 
the electric held in a horizontal cross-section through 
the grating region at z = h/2 is depicted in hgure 
As in the ID case, these maps show that the held 
exhibits spatial oscillations with smaller amplitude as 
compared to the variations of the held especially 
at ^-aligned interfaces (outlined with blue dashed lines 
in hgure 1^. 

5. Out-of-plane normal vector fields for 
oblique diffraction gratings 

Oblique dihraction gratings as those shown in hg- 
ures Sb) and [^d) are modeled in the RCWA within 
the staircase approximation along a direction perpen- 
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Figure 9. a) Spatial distribution of the dominant component 
of the electric field, \Ex\^ in the 2D grating region at 2 ; = h/2, 
determined by using the conventional RCWA and N = 27. b) 
Spatial distribution of the dominant component of the electric 
near-field, \Ex\^ determined for the same grating parameters as 
in a) but using the improved algorithm. The blue vertical line 
was added for clarity and merely separates the two plots. 


dicular onto the grating plane. Each computational 
slice is assumed to be a 2 ;-independent structure, for 
which the eigenmodes can be found numerically by 
solving equation 0 . The field in each computational 
layer is then found by using the boundary conditions at 
the top and bottom interfaces of the grating, employ¬ 
ing a numerically stable 5-matrix formulation. The 
validity of this staircase approximation for ID-periodic 
gratings under TE-polarization has been proven in . 
In the context of the NVE formulation of RCWA for 
oblique ID-periodic structures, it has been found that 
the use of an out-of-plane NVE, i.e. N = 0, Nz)^^ 

is beneficial for the far-held convergence speed [^[M] . 
In this section we will study the relation between the 
accurate held formulation ( 10 ) and the near-held con¬ 
vergence speed for slanted ID-periodic structures. 

The situation for oblique 2D-periodic structures 
has not yet been explored in the context of RCWA. 
However, based on the results presented so far, one 
can conjecture that both the near- and far-held 
convergence of RCWA can be improved by using an 
out-of-plane 3D NVE, N = Ny^ 7 ^ 0, and 
that the near-held evaluation would be more accurate 
as well. The validity of this supposition will be 
explored in the second part of this section. 


5.1. Analysis of slanted ID-periodic binary diffraction 
gratings 

In order to analyze oblique ID-periodic structures, 
we consider the slanted binary grating depicted in 
hgure Bb). The period of the grating is A = 1 pm, 
the hlling factor, p = 0.5, the height, h = 0.25 pm, 
and the slanting angle is 0 = 7 r/ 4 . Only the gold 
grating is considered in this section as this would be 
the most challenging case. If the unit cell is assumed 




Figure 10. Computational results for the slanted ID binary 
grating shown in figure [^b). a) Far-field self-error ep'{N, M) vs. 
N and M corresponding to the in-plane NVF formulation, b) 
Far-field self-error eF(N^M) corresponding to the out-of-plane 
NVF formulation, c) Near-field self convergence 
corresponding to in-plane NVF. d) Near-field self convergence 
corresponding to out-of-plane NVF. 


to extend from x = —A/2 to x = A/2 and the center 
of the binary grating is set to be x = 0, z = h/2, 
a suitable out-of-plane NVE is given by N(x,^, 2 :) = 
sign {x — z) \/2 (1, 0,1)^. 

The two formulations of the RCWA compared in 
this section are the in-plane NVE, which is used in 
conventional RCWA together with the conventional 
field evaluation (|^, and the out-of-plane NVE, 
N, combined with the improved field evaluation 
formulation (10). In contrast to the results presented 
in the previous sections, the two formulations yield 
different results for both the near- and far-held 
quantities. Moreover, for the sake of the clarity of the 
presentation, all physical quantities corresponding to 
the out-of-plane NVE formulation are denoted with a 
tilde symbol. 

Numerical results for increasing number of har¬ 
monics, V = 2,..., 320, and number of computational 
layers, M = 2,..., 256, are presented in hgure It 
can be inferred from this hgure that the in-plane formu¬ 
lation requires both a high number of ES coefficients, 
V, and layers, M, to achieve convergence to a result of 
ppref _ 0.28788, whereas the out-of-plane formulation 
yields fast convergence to = i^(320,256) _ 0.28837 
with respect to V, as per hgnres [l^a) and p^b), re¬ 
spectively. This behavior is in agreement with the hnd- 
ings reported in 24 . The convergence of the calculated 
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Figure 11. a) Spatial distribution of the electric near-field, 
1^^40,256) I ^ determined using the conventional in-plane NVF. 
b) Spatial distribution of the electric near-field, 
determined for the same grating parameters as in a) but using 
the modified out-of-plane NVF formulation. 




near-field, illustrated in figures [1 Q]^c) and [Tol^d), ex¬ 
hibits similar features. Specifically, the in-plane NVF 
formulation requires both high N and M to achieve a 
small self-error of } = 4.7 • 10“^, whereas 

this self-error can already be achieved with N = 10, 
M = 256 in the out-of-plane formulation. This clearly 
demonstrates a drastically improved efficiency to the 
calculation of the near-held of oblique diffraction grat¬ 
ings of the approach based on the combination of out- 
of-plane NVF and the accurate near-held formulation. 
The highly improved near-held prohle is illustrated in 
hgnrepTj^b), which exhibits no unphysical oscillations 
near the gold-vacuum interface. This is in sharp con¬ 
trast to the conventional held evaluation of the in-plane 
formulation (cf. hgure [lQ]^a)), which clearly suhers 
from spurious oscillations. 

5.2. Analysis of slanted 2D-periodic cylindrical 
diffraction gratings 

In this section we investigate the efficiency of using 
the accurate held evaluation and the out-of-plane 
NVF formulation to model a challenging, slanted 2D- 
periodic dihraction grating. The grating with periods, 
Ai = A 2 = A = 1 pm, is schematically depicted 
in hgure Bd) and consists of a cylindrical rod with 
radius r = 0.3 pm and height h = 0.125 pm, which is 
slanted by 0 = 7r/4 along the x-axis. Again, only the 
gold grating is considered in this section. Finally, the 
incident plane wave is impinging onto the grating along 
the normal direction, is polarized along the x-axis, and 
has a wavelength of A = 2 pm. 

The rehection coefficient calculated for N = 
3,..., 19 harmonics and M = 2,..., 32 layers is shown 


Figure 12. Computational results for the slanted ID binary 
grating shown in figure [^d). a) Refiection coefficient R vs. N 
and M, determined using the in-plane NVF formulation, b) 
Refiection coefficient R vs. N and M, determined using the 
out-of-plane NVF formulation, c), d) Near-field convergence of 
conventional and modified formulation, respectively. 

in hgure [T^a) and hgure [T^b) and was determined 
by using the conventional in-plane NVF and the 
out-of-plane NVF formulation, respectively. It can 
be seen that both approaches converge rather slow, 
neither one achieve convergence even for the highest 
considered values of M = 32 and V = 19. Moreover, 
in order to characterize the error of the near-held 
calculations in the two formulations, the self-error 
with respect to the reference solutions obtained with 
V = 19 and M = 32 is presented in hgures [T^c) 
an d@d). The in-plane formulation achieves a relative 
self-error of = 0.625, whereas for the same 

values of N and M, the modihed held evaluation 
in conjunction with the out-of-plane NVF achieves a 
substantially lower self-error of = 0.261. 

It has to be stressed, that the necessary accuracy 
for full convergence could not be achieved in our 
simulations. The evolution of the computational 
results for the shown values of V < 19 and M < 32 
however can be interpreted in favor of the out-of-plane 
NVF formulation, due to the lower near-held self¬ 
error }. It can be supposed 

that future simulations with hner discretization will 
reveal the practical beneht of the out-of-plane NVF 
formulation in conjunction with the modihed held 
formulation for 2D-periodic slanted structures, similar 
to the case of ID-periodic slanted dihraction gratings. 
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6. Conclusions 

To summarize, we have analyzed the numerical 
near-held calculated using the RCWA and identihed 
the Gibb’s phenomenon as the main reason for 
their slow convergence and the spurious oscillations 
from which they suffer. As a solution to these 
dehciencies of RCWA, we proposed an improvement 
of this methods that can be applied for modeling 
arbitrary diffraction gratings and, more generally, 
periodic optical structures. The modihed formulation 
signihcantly improves the accuracy of ID-RCWA 
calculations, for both straight and slanted gratings, 
where it speeds up convergence and removes the 
numerical artifacts from the calculated near-helds. The 
accuracy of 2D-periodic grating simulations can be 
enhanced, however, to a lesser extent than in the ID 
case. The reduced performance in 2D can be attributed 
to the discontinuity and non-exactness of the numerical 
NVF, which is at the core of the modified formulation. 
Therefore, it might be fruitful to investigate more 
elaborate NVF-formulations and their suitability for 
near-held calculations, such as a complex valued NVF 
[ 41 ] , which unlike the NVF used here is continuous 
everywhere in the grating region. 

We expect that the proposed modihcation of the 
RCWA method will greatly advance its computational 
capabilities, especially for ID periodic optical struc¬ 
tures. In particular, this improved method could 
prove instrumental to accurate modeling of periodic 
plasmonic structures, diffraction gratings, and surface- 
nonlinear devices, namely to simulation of physical sys¬ 
tems whose functionality rely on the electromagnetic 
near-held at interfaces. 
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